#' ---
#' title: "More donors, more democracy: Appendix B: The instrumental variable approach"
#' subtitle: "Placebo test for spurious instruments"
#' output:
#'   pdf_document:
#'     toc: true
#'     number_sections: true
#' author: Sebastian Ziaja
#' date: "This version: 3 November 2018"
#' ---

#+ setup

numboot <- 10 # set to 1000
startYear <- 1994
endYear <- 2013

workdir <- ""
# workdir <- ifelse(basename(getwd()) == "aiddem_new", "", "../")
source(paste0(workdir, "scripts/load_packages.R"))
source(paste0(workdir, "scripts//extract.custom.R"))
source(paste0(workdir, "scripts//functions.R"))
source(paste0(workdir, "scripts//formulas.R"))
source(paste0(workdir, "scripts//variable_labels.R"))
source(paste0(workdir, "scripts//getData.R"))
source(paste0(workdir, "scripts//aggData.R"))

# determine variable set for core sample generation
vars_core_set <- all.vars(as.formula(bf[[4]]))

for(ctrlvars in c(#T,   # include control variables?
                  F)) {
for(AB in c(#"XX", "AX",  # randmize which part of dyad?
            "XB")) {
for(iax in c(#T,  # use interaction specification?
             F)) {
    for(wax in c(1 # yearly data or several-year averages?
                 # 2, 3
                 )) {
        w <- wax
        for(iam in c(#T, # include multilateral donors?
                     F)) {
            for(rax in c(#T,
                         F)) { # restrict sample to non-electoral regimes

# Dataset parameters (set via loop, do not modify):
includeMultilateralDonors <- iam
interactX <- iax
wx <- wax

                cat("SETTINGS:  ", wx, "_", numboot,
                  ifelse(ctrlvars, "_ctrl", ""),
                  ifelse(includeMultilateralDonors, "_m", ""),
                  ifelse(interactX, "_IA", ""),
                  ifelse(rax, "_reg0", ""), AB, "\n\n")

# dx <- getData(useCache = FALSE, useCacheCov = FALSE, advx = "3.1")

# d1 <- aggData(w = 1)
# d2 <- aggData(w = 2)
# d3 <- aggData(w = 3)
# d4 <- aggData(w = 4)


dxp <- getData(useCache = FALSE, useCacheCov = FALSE, advx = "3.1",
               placebo = TRUE, bootx = numboot)
dXp <- aggData(dat = dxp, w = wx,
               minobs = ifelse(w == 1, 10,
                               ifelse(w == 2, 5, 3))
)

if(rax) {
    dp <- dXp %>% filter(smplR == 1)
} else {
    dp <- dXp %>% filter(smpl == 1)
}

if(ctrlvars) { bx <- bf[[3]] } else {
bx <- sub("l.pop.log.r + l.gdpcap.log.r + l.war25", "1", bf[[3]], fixed = TRUE)
}
if(interactX) {
    bx <- sub("l.CMgnh", "IA_l", bx)
}
om1 <- felm(as.formula(bx), data = dp)
summary(om1)
ncoef <- length(om1$coef) - (length(om1$endovars) - 1)
oest <- om1$coeff[[ncoef]]

n <- numboot
r <- vector("numeric", n)
for(i in 1:n) {
    r[i] <- felm(as.formula(
        sub("l.ZwvCMgwh94", paste0("Zg", i, "_l"), bx)
    ), data = dp)$coeff[[ncoef]]
}
dens <- density(r, n = 10^5)

# plot(dens$x, dens$y, type = "l")
# abline(v = 0, lty = 2)
# abline(v = oest, lwd = 2, col = "green")
# abline(v = mean(r), lwd = 2, col = "red", lty = 2)
# abline(v = median(r), lwd = 2, col = "red", lty = 3)

{
pdf(file = paste0(workdir, "graphs/randomplacebo_w", wx, "_", numboot,
                  ifelse(ctrlvars, "_ctrl", ""),
                  ifelse(includeMultilateralDonors, "_m", ""),
                  ifelse(interactX, "_IA", ""),
                  ifelse(rax, "_reg0", ""), AB,
                  ".pdf"),
    width = 8, height = 5)
par(mfrow = c(1, 1))
plot(dens$x, dens$y, type = "l", xlim = if(interactX) { c(-.3, .7) } else { c(-1, 2.5) }
     # , xaxt = "n"
     , xlab = "IV estimate from baseline model"
     , ylab = "density"
     )
abline(v = 0, lty = 2)
abline(v = oest, lwd = 2, col = "red")
abline(v = mean(r), lwd = 2, col = "red", lty = 2)
abline(v = median(r), lwd = 2, col = "red", lty = 3)

pv <- quantile(r, c(.05, .1, .9, .95))

# polygon(c(dens$x[dens$x <= pv[2]], pv[2]),
#         c(dens$y[dens$x <= pv[2]], 0), col = "gray")
# polygon(c(dens$x[dens$x >= pv[3]], pv[3]),
#         c(dens$y[dens$x >= pv[3]], 0), col = "gray")
polygon(c(dens$x[dens$x <= pv[1]], pv[1]),
        c(dens$y[dens$x <= pv[1]], 0), col = "gray")
polygon(c(dens$x[dens$x >= pv[4]], pv[4]),
        c(dens$y[dens$x >= pv[4]], 0), col = "gray")

legend("topleft",
       title = "Estimates",
       legend = c(paste0("Original (",
                         prettyNum(oest, digits = 2), ")"),
                  paste0("Mean placebo (",
                         prettyNum(mean(r), digits = 2), ")"),
                  paste0("Median placebo (",
                         prettyNum(median(r), digits = 2), ")")),
       col = "red", cex = .75,
       lty = 1:3, lwd = 2
       )

legend("topright",
       title = "Bottom and top draws",
       legend = c(
           paste0("5 percent [",
                         prettyNum(quantile(r, .05), digits = 2), "; ",
                         prettyNum(quantile(r, .95), digits = 2),
                         "]")
                  # , paste0("10 percent [",
                  #        prettyNum(quantile(r, .1), digits = 2), "; ",
                  #        prettyNum(quantile(r, .9), digits = 2),
                  #        "]")
                  ),
       cex = .75,
       fill = c(#"black",
           "gray")
       )
dev.off()
}

}}}}}}
